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PREFACE 


The  goal  of  the  work  reported  here  for  the  Underground  Technology  Program  (UTP)-  Three- 
Dimensional  Calculations,  is  to  predict  the  dynamic  response  and  damage  to  reinforced  tun¬ 
nels  in  intact  and  jointed  rock  subjected  to  ground  shock  from  nuclear  weapons.  The  finite 
element  calculations  were  performed  for  the  Defense  Nuclear  Agency  (dna)  by  Ms.  Yvonne 
Murray  of  APTEK  Inc.,  and  consultant  Dr.  Leonard  Schwer  of  Schwer  Engineering  and  Con¬ 
sulting  Services  (SE&CS)  under  contract  DNA001-91-C-0140.  The  APTEK  Program  Manager 
was  Dr.  Eugene  Fitzgerald.  The  DNA  technical  monitor  was  Dr.  Paul  Senseny. 

This  report  is  a  concise  summary  of  APTEK ’s  effort,  with  emphasis  on  the  calculations 
performed  by  Ms.  Murray.  These  calculations  include  dynamic  tunnel  test  and  stack  of  brick 
simulations,  with  and  without  modeling  rate  effects  and  damage.  This  report  highlights  what 
we  learned  and  provides  recommendations  for  future  efforts.  More  thorough  summaries  of 
our  intact  rock  and  joint  models  are  documented  in  an  APTEK  technical  report  (Murray  and 
Lindberg,  1996). 

Lessons  learned  from  Dr.  Schwer’s  simulations  are  documented  by  Thacker,  Riha,  and  Schwer 
(1995).  These  calculations  include  static  and  preliminary  dynamic  tunnel  test  simulations. 
Emphasis  is  on  rate-dependent  and  rate-independent  simulations,  without  modeling  damage. 
This  document  also  compares  numerical  simulations  using  the  Sandler  and  Rubin  (1979)  cap 
model  with  the  smooth-cap  model  originally  developed  by  Pelesonne  (1989)  and  enhanced 
by  APTEK.  Both  Ms.  Murray’s  and  Dr.  Schwer’s  calculations  were  performed  as  part  of  the 
Precision  Test  Modeler’s  (PTM)  benchmarking  and  validation  effort. 

In  addition  to  the  modeling  and  simulation  effort  described  in  this  report,  tunnel  kill  algo¬ 
rithms  for  the  Lethality  Analysis  of  Buried  Structures  (laBS)  software  package  were  devel¬ 
oped  by  consultant  Dr.  Herbert  Lindberg  of  LCE  Software/Engineering.  A  technical  reference 
which  gives  details  of  the  closed-form  tunnel  kill  algorithms  and  their  technical  justification 
is  currently  available  through  APTEK  (Murray  and  Lindberg,  1996),  and  will  be  published 
as  a  DNA  topical  report  at  a  later  date. 
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CONVERSION  TABLE 


Conversion  factors  for  U.S.  Customary  to  metric  (SI)  units  of  measurement. 
MULTIPLY  .  BY  ^  —  TO  GET 

TO  GET  < —  < —  BY  < —  < —  DIVIDE 


angstrom 

1.000  000  X  E-10 

meters  (m) 

atmosphere  (normal) 

1.013  25  X  E  +2 

kilo  pascal  (kPa) 

bar 

1.000  000  X  E  +2 

kilo  pascal  (kPa) 

barn 

1.000  000  X  E  -28 

meter^  (m^) 

British  thermal  unit  (thermochemical) 

1.054  350  X  E  -1-3 

joule  (J) 

calorie  (thermochemical) 

4.184  000 

joule  (J) 

cal  (thermochemical)/cm^ 

4.184  000  X  E  -2 

megajoule/m^  (MJ/m^) 

curie 

3.700  000  X  E  -hi 

giga  becquerel  (GBq)* 

degree  (angle) 

1.745  329  X  E-2 

radian  (rad) 

degree  Fahrenheit 

r^  =  (t°/  +  459.67)/1.8 

degree  kelvin  (K) 

electron  volt 

1.602  19  X  E-19 

joule  (J) 

erg 

1.000  000  X  E-7 

joule  (J) 

erg/second 

1.000  000  X  E-7 

watt  (W) 

foot 

3.048  000  X  E  -1 

meter  (m) 

foot-pound-force 

1.355  818 

joule  (J) 

gallon  (U.S.  liquid) 

3.785  412  X  E  -3 

meter^  (m^) 

1  inch 

2.540  000  X  E  -2 

meter  (m) 

jerk 

1.000  000  X  E  -4-9 

joule  (J) 

joule/kilogram  (J/kg) 

1.000  000 

Gray  (Gy)  ** 

(radiation  dose  absorbed) 

kilotons 

4.183 

terajoules 

kip  (1000  Ibf) 

4.448  222  x  E  -1-3 

newton  (N) 

kip/inch“  (ksi) 

6.894  757  x  E  -t-3 

kilo  pascal  (kPa) 

ktap 

1.000  000  X  E  -1-2 

newton-  second/m^  (N-s/m^) 

micron 

1.000  000  X  E  -6 

meter  (m) 

mil 

2.540  000  X  E  -5 

meter  (m) 

mile  (international) 

1.609  344  X  E  -1-3 

meter  (m) 

ounce 

2.834  952  x  E  -2 

kilogram  (kg) 

pound-force  (Ibf  avoirdupois) 

4.448  222 

newton  (N) 

pound-force  inch 

1.129  848  X  E  -1 

newton-meter  (N*m) 

pound-force/inch 

1.7.51  268  X  E  -1-2 

newton/meter  (N/m) 

pound-force/foot“ 

4.788  026  X  E  -2 

kilo  pascal  (kPa) 

pound-force/inch“  (psi) 

6.894  757 

kilo  pascal  (kPa) 

pound-mass  (Ibm  avoirdupois) 

4  . .535  924  x  E  -1 

kilogram  (kg) 

pound-mass-foot“  (moment  of  inertia) 

4.214  011  X  E  -2 

kilogram-meter^  (kg#m^) 

pound-mass/foot^ 

1.601  846  X  E  -bl 

kilogram/meter^  (kg/m^) 

rad  (radiation  dose  absorbed) 

1.000  000  X  E-2 

Gray  (Gy)** 

roentgen 

2.579  760  X  E  -4 

coulomb/kilogram  (C/kg) 

shake 

1.000  000  X  E -8 

second  (s) 

slug 

1.459  390  X  E  -El 

kilogram  (kg) 

torr  (mm  Hg,  O^C) 

1.333  22  X  E-1 

kilo  pascal  (kPa) 

*  The  becquerel  (Bq)  is  the  SI  unit  of  radioactivity;  1  Bq  =  1  event/s. 
♦♦The  Gray  (Gy)  is  the  SI  unit  of  absorbed  radiation. 
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SECTION  1 

INTACT  ROCK  MODEL  AND  SWAT  SIMULATIONS 


1.1  INTACT  ROCK  MODEL  OVERVIEW. 

Performing  calculations  with  a  first  principles  code  requires  implementation  of  a  comprehen¬ 
sive  material  model  for  rock  and  other  geological  materials.  We  implemented  a  smooth-cap 
model  with  damage  into  the  nonlinear  finite  element  code  DYNA3D  (Whirley,  1993)  to  an¬ 
alyze  the  static  and  dynamic  response  of  tunnels  in  intact  rock.  The  original  basis  of  our 
geological  material  model  is  a  two-invariant,  smooth-cap  elasto-plastic  model  (Pelesonne, 
1989).  We  added  numerous  features  to  this  model  to  improve  fits  to  standard  laboratory 
test  data.  These  improvements  include  a  scalar  anisotropic  damage  formulation  for  model¬ 
ing  strain-softening  and  modulus  reduction,  viscoplastic/viscodamage  formulations  and  rate- 
shifted  damage  thresholds  for  modeling  high  strain  rate  effects,  an  element  length  scale  to  reg¬ 
ularize  the  damage  dissipation  to  be  mesh  size  independent,  and  a  three-invariant  plasticity^ 
formulation  to  simultaneously  fit  triaxial  tension  and  compression  data.  The  model  captures 
the  essential  features  of  geological  material  behavior;  shear  enhanced  compaction,  dilatency, 
pre-peak  hardening,  post-peak  softening,  modulus  reduction,  and  localized  damage  accu¬ 
mulation.  Good  fits  of  the  model  to  quasi-static  and  high  strain  rate  data  are  given  in 
Appendix  A. 

1.2  SWAT  TEST  OVERVIEW. 

SRI  conducted  spherical  wave  tests  (s\v.\T)  on  aluminum  lined  tunnels  in  limestone  cylinders 
(Klopp  et.  ai  1995),  as  shown  in  Figure  1-1.  The  objective  of  these  tests  was  to  increase 
our  understanding  of  dynamic  tunnel  closuix'  mechanisms  and  to  provide  data  for  model 
verification.  Numerous  tunnels  were  simultaiu'ously  tested  and  located  at  various  distances 
the  center  of  the  PETN  charge.  SRI  n'corded  free-field  velocity  and  stress  histories,  and  the 
crown-invert  and  springline  closures. 

Using  our  physically-based  model.  w<‘  (•alculat('d  tunnel  closures,  free-held  stress  histories, 
and  free-held  velocity  histories  for  comparison  with  the  recorded  SWAT  data.  Our  hnite 
element  mesh  with  two  tunnels  modeled  at  1  1.5  cm  and  19.2  cm  is  shown  in  Figure  1-2. 
Instead  of  modeling  the  explosive  charne.  we  ap|)lied  the  measured  velocity  history  shown 
in  Figure  1-3  to  the  inside  of  the  7.5  cm  radius  cavity.  This  is  the  velocity  recorded  by  gage 
PVl  at  that  location. 

SRI  measured  the  closures  on  the  liner,  not  the  tunnel.  Hence  we  report  predicted  liner 
closures,  rather  than  tunnel  closures.  We  modeled  a  frictionless  sliding  interface  between 
the  tunnel  and  liner  which  allows  tlu'  liner  to  separate  from  the  tunnel.  The  0.79-mm  thick 
liner  is  made  of  annealed  aluminum  and  modeled  as  an  elastic-plastic  material.  SRI  did  not 
report  whether  the  liner  separated  from  the  tunnel  during  the  test. 

^This  formulation  was  iinplementecl.  in  part,  on  d.na’.s  Conventional  Weapons  Effects  Program  and  is 
reported  by  Murray  and  Lewis  (1995)  and  .Scliwer  and  Murray  (1994). 
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Figure  1-1.  SWAT  dynamic  test  configuration. 

Numerous  calculations  were  performed  and  compared  with  test  data  to  determine  the  effects 
of  modeling  strain-rate  dependence  and  damage.  Comparisons  include  liner  closures  at  the 
14.5  cm  and  19.2  cm  tunnel  locations,  free  field  velocity  histories  at  8.33  cm  (PV2),  13.33  cm 
(PV5).  18.33  cm  (PV8).  and  23.33  cm  (PVll)  from  the  center  of  the  charge,  and  free  field 
stress  histories  at  14.5  cm  and  19.2  cm. 

1.3  SWAT  C0MPUT.4TI0NAL  SUMMARY. 

Overall,  the  rate-independent  calculation  without  damage  simulates  the  free-field  velocity 
and  stress  histories  reasonably  well,  although  the  calculated  peak  velocities  and  late-time 
stress  histories  are  slightly  lowered  than  measured.  This  is  demonstrated  in  Figure  1-4. 
However,  the  calculated  crown-invert  closure  at  14.5  cm  is  about  25%  lower  than  measured^. 
In  addition,  the  calculated  crown-invert  closure  at  19.2  cm  is  about  50%  lower  than  measured, 
and  does  not  exhibit  the  rebound  of  the  mecisured  closure. 

One  possible  explanation  for  the  difference  between  the  measured  and  calculated  crown 
invert  closures  at  19.2  cm  is  that  the  liner  separated  from  the  tunnel  during  the  test,  then 

-In  Figures  l-4a  and  l-4b  we  show  both  the  crown-invert  and  springline  closures  on  a  single  plot.  Using 
the  tunnel  located  at  14.5  cm  as  an  example,  the  calculated  crown-invert  closure  is  about  5%  while  the 
measured  closure  is  about  6.791.  The  calculated  springline  closure  is  about  1%  while  the  measured  closure 
is  about  2%. 
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Figure  1-2.  Close-up  of  the  SWAT-2  finite  element  mesh  near  the  source  and  tunnels. 


Figure  1-3.  Measured  SWAT  velocity  applied  to  inner  cavity  in  finite  element  analysis. 
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(a)  Closures  at  the  14.5  cm  tunnel. 
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(b)  Closures  at  the  19.2  cm  tunnel. 
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(c)  Velocity  histories  at  four  gage  locations. 
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(d)  Stress  histories  at  two  gage  locations. 


Figure  1-4.  The  rate-independent  calculation  without  damage  underpredicts  the  measured 
crown-invert  closures. 
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rebounded  to  join  the  tunnel  at  late  time  (around  0  25  msec).  The  li;-;  ;r  did  not  separate  from 
the  tunnel  during  the  calculation,  although  the  r.-alcuiafced  tunnel  clci  ure  is  in  agreement  with 
one  late  time  measurement.  Some  preliminary  calculations  with  a  less  refined  mesh  indicated 
separation  at  this  tunnel  location. 

Adding  viscoplasticity  (rate-dependence)  to  the  calculation  without  damage  improves  com¬ 
parisons  between  the  measured  and  calculated  velocity  histories,  and  between  the  measured 
and  calculated  crown-invert  closure  at  the  14.5  tunnel  location.  This  is  demonstrated  in  Fig¬ 
ure  1-5.  However,  viscoplasticity  has  little  effect  on  the  crown-invert  closure  calculated  at 
the  19.2  cm  tunnel.  Modeling  viscoplasticity  also  opens  the  springline  closures  at  both  tun¬ 
nel  locations,  in  lack  of  agreement  with  the  measured  histories.  We  note  that  the  measured 
springlines  open  slightly  at  early  time  (around  0.5-0. 6  msec)  before  closing.  Our  calculation 
starts  out  with  the  right  springline  opening,  but  continues  to  open  instead  of  close.  This 
suggests  that  some  mechanism,  not  currently  modeled,  initiates  the  closing.  However,  SRI 
reports  that  the  resolution  of  the  gages  is  0.3%  of  the  tunnel  closure  (Klopp  et.  al,  1995), 
and  the  initial  opening/closing  is  on  the  order  of  this  resolution. 

Modeling  viscoplasticity  near  the  source,  but  not  near  the  tunnels,  provides  excellent  agree¬ 
ment  between  the  measured  and  calculated  histories  at  almost  all  locations.  This  is  demon¬ 
strated  in  Figure  1-6.  One  exception  is  the  crown-invert  closure  at  the  19.2  cm  tunnel  in 
which  the  calculated  closure  agrees  with  the  measured  closure  at  late  time,  but  not  through¬ 
out  the  entire  history.  The  other  exception  is  the  early-time  opening/closing  behavior. 

It  is  interesting  to  note  the  free-field  velocity  and  stress  histories  calculated  in  Figures  1-5 
and  1-6  are  only  subtlety  different,  yet  the  calculated  springline  closures  are  significantly 
different  (of  opposite  sign).  By  not  modeling  of  rate  effects  near  the  tunnel,  we  were  able  to 
change  the  calculated  springline  deformation  from  opening  to  closing  at  the  14.2  cm  location. 

To  help  understand  these  correlations,  it  is  necessary  to  examine  the  implementation  of  the 
viscoplastic  model  and  how  it  was  fit  to  high  strain  rate  data.  The  behavior  of  the  viscoplastic 
model  and  fits  to  high  strain  rate  data  are  thoroughly  discussed  in  the  Section  1.4,  so  only 
the  pertinent  issues  are  discussed  here. 

Effectively,  we  model  different  rate  effects  on  the  cap  than  on  the  shear  failure  surface. 
This  is  done  by  fitting  our  viscoplastic  model  to  two  sets  of  data,  with  a  smooth  transition 
between  fits.  These  data  are  high  pressure  data  (on  the  cap)  and  uniaxial  stress  data  (on 
shear  failure  surface).  Examples  of  high  pressure  data  are  Hugcniot  data  or  SWAT  free-field 
data.  We  set  the  cap  (pressure-dependent)  viscoplastic  parameters  to  get  good  agreement 
with  the  SWAT  PV  2  velocity  gage  history^.  We  approximately  set  the  the  shear  failure  surface 
parameters  to  fit  Hopkinson  split  pressure  bar  data  (compressive).  Based  on  a  general  review 
of  Hopkinson  bar  data  available  in  the  literature,  including  that  lor  Solenhofen  limestone 
(Salem  limestone  data  was  not  available),  we  assumed  that  strength  enhancement  increases 

®We  also  fit  the  pressure-dependent  viscoplastic  parameters  to  Salem  Limestone  Hugoniot  data  available 
in  the  literature,  as  discussed  in  Appendix  A.  However,  the  swat  velocity  his'ory  correlations,  using  the 
Hugoniot  fit,  were  not  as  good  as  those  shown  in  Figure  1-5,  using  the  pv2  gag-  fit,  so  they  are  not  shown 
here.  Similarly,  modeling  viscoplasticity  near  the  tunnel  (with  the  Hugoniot  fit)  oj.-ened  the  springline  slightly, 
consistent  with  the  trend  already  noted. 


5 


Velocity  (m/s)  closures  (%) 


14.5  cm 


Time  (usee) 

(a)  Closures  at  the  14.5  cm  tunnel. 
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(b)  Closures  at  the  19.2  cm  tunnel. 
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(d)  Stress  histories  at  two  gage  locations. 


Figure  1-5.  The  rate-dependent  calculation  without  damage  is  in  good  agreement  with  the 
measured  velocity  and  stress  histories,  but  predicts  springline  opening  instead 
of  closing. 
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(a)  Closures  at  the  14.5  cm  tunnel. 
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(b)  Closures  at  the  19.2  cm  tunnel. 
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(c)  Velocity  histories  at  four  gage  locations.  (d)  Stress  histories  at  two  gage  locations. 


Figure  1-6.  Modeling  rate  effects  near  the  source,  but  not  near  the  tunnels  provides  good 
agreement  with  most  of  the  measured  velocit5\  stress,  and  closure  histories. 


by  a  factor  of  two  at  high  strain  rates  near  10'*.  Hopkinson  bar  data  for  Solenhofen  limestone 
is  given  in  Figure  1-12  of  Section  1.4. 

The  observation  that  we  calculate  the  free-field  velocity  and  stress  histories  adequately  in 
Figure  1-5,  but  not  the  springline  closures,  suggests  that  we  model  high  pressure  (cap)  rate 
effects  adequately,  but  not  low  pressure  (shear  failure  surface)  rate  effects.  Near  the  source, 
the  pressure  is  high  and  the  stress  state  lies  on  the  cap.  Near  the  tunnels,  the  pressure  is 
low  and  the  stress  path  encounters  the  shear  failure  surface.  We  attempted  one  calculation 
in  which  we  modeled  minimal  rate  effects  on  the  shear  failure  surface;  we  assumed  that 
strength  enhancement  increase  by  10%  at  high  strain  rates  near  10^,  but  this  calculation  did 
not  improve  correlations  with  the  springline  measurement.  The  only  way  we  achieved  good 
correlations  with  the  springline  measurement  was  to  turn  viscoplasticity  off  completely  near 
the  tunnels,  as  previously  demonstrated  in  Figure  1-6. 

Two  alternative,  but  tentative  conclusions  could  be  drawn  here.  The  first  conclusion  is  that 
our  fits^  and/or  viscoplastic  formulation  need  tuning.  This  conclusion  is  partly  supported  by 
the  SWAT  correlations  of  other  PTM  group  calculators,  specifically  Weidlinger  &:  Associates 
(wa)  and  Titan  Research  Corporation  (trt).  By  adjusting  their  rate- dependent  fits  on  the 
cap,  and  modeling  minimal  rate  effects  on  the  shear  failure  surface,  both  calculators  were 
able  to  simulate  springline  closure  at  the  14.2  cm  tunnel.  However,  at  the  time  this  report 
was  written,  no  PTM  calculator  was  able  to  simulate  the  the  crown-invert  rebound  and  slight 
springline  closure  measured  at  the  19.2  cm  tunnel.  In  addition,  no  calculator  demonstrated 
that  tuning  the  viscoplastic  model  to  the  SWAT  data  provides  a  good  predictive  capability 
for  other  tunnel  tests,  because  no  such  predictions  were  made.  Future  efforts  should  include 
additional  tunnel  response  predictions,  particularly  for  full  scale  field  tests. 

The  second  possible  conclusion  is  that  we  need  to  examine  alternatives  to  viscoplasticity,  such 
as  rate-shifted  plasticity  surfaces  or  better  dynamic  stress-strain  relations.  When  generating 
high  strain  rate  Hugoniot  data  like  that  shown  in  Figure  A-3a  of  Appendix  A,  peak  stress 
and  strain  values  are  derived  from  measured  shock  and  particle  velocity  histories  (Larson 
and  Anderson,  1979).  This  derivation  is  usually  based  conservation  laws  for  uniaxial  flow, 
or  the  Rankine-Hugoniot  jump  equations.  Viscoplasticity  is  not  assumed  in  the  derived 
stress/strain  values,  although  pore  collapse  rate  effects,  if  any,  may  affect  the  measured 
particle  velocity  histories. 

All  PTM  calculators,  including  APTEK,  used  a  constant  bulk  modulus  model  for  the  pressure- 
volumetric  strain  relation,  as  follows: 


P  =  Kie,-el)  (1.1) 

where  P  is  the  pressure,  K  is  the  bulk  modulus,  is  the  total  volumetric  strain,  and  ej  is 
the  plastic  volumetric  strain.  Once  all  pores  are  collapsed,  ej  remains  constant  its  maximum 
value,  so  Equation  1.1  simulates  a  /in ear  relation  between  the  pressure  and  volumetric  strain, 
at  least  at  high  pressure. 

■^One  difference  between  our  fits  to  laboratory  data,  and  those  of  the  other  PTM  calculators,  is  that  our 
cap  was  initially  located  at  a  lower  pressure  than  those  of  the  other  calculators.  This  could  make  a  difference 
when  modeling  pressure-dependent  rate  eflfects. 
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(a)  High  stress  behavior.  (b)  Low  stress  behavior. 


Figure  1-7.  The  Mie-Gruenisen  model  predicts  the  stiffening  behavior  of  geological  materials 
at  high  stress  better  than  the  constant  bulk  modulus  model. 

For  consistency  in  interpreting  and  fitting  Hugoniot  data,  we  suggest  evaluating  an  alter¬ 
native  pressure- volumetric  strain  relation,  based  on  the  Rankine-Hugoniot  relations  for  a 
travelling  shock  wave,  and  the  assumption  of  a  linear  shock- velocity  versus  particle  velocity 
relationship,  as  follows: 

^  pc^(t>, 

{1-S4>cY 

where  p  is  the  density,  c  is  the  bulk  sound  speed.  5  is  the  slope  of  the  linear  fit  to  the  shock 
velocity  versus  particle  velocity  data,  and  Oc  is  the  relative  change  in  volume  of  the  solid 
particles.  0c  depends  on  the  total  and  plastic  volume  strains,  as  discussed  by  Murray  and 
Lewis  (1995).  This  relation  was  incorporated  into  an  explicit  pore  collapse  model  with  a 
Mie-Gruniesen  equation  of  state  on  DNA  s  CVVE  program.  The  three  user-supplied  material 
parameters  for  the  Gruenisen  model  are  c  (or  K  =  pc^)  and  S  and  the  Gruenisen  parameter 
r,  compared  with  K  for  the  constant  modulus  model. 

Although  the  Mie-Gruniesen  model  was  not  used  on  the  current  UTP  program,  a  fit  to  soil 
data  was  previously  made  on  the  CWE  program  which  illustrates  the  general  behavior  of  the 
model.  The  low  stress  behaviors  of  the  constant  modulus  model  in  Equation  1.1  and  the  Mie- 
Gruniesen  model  in  Equation  1.2  are  nearly  identical,  but  the  high  stress  behaviors  differ,  as 
shown  in  Figure  1-7.  Equation  1.2  simulates  a  nonlinear  relation  between  the  pressure  and 
volumetric  strain,  which  stiffens  at  high  stress  levels. 

Future  efforts  should  examine  the  specific  behavior  of  the  Mie-Gruniesen  and  constant  bulk 
modulus  models  for  Salem  limestone,  using  approximate  values  for  S  and  F  available  for  lime¬ 
stone  (or  rock)  in  the  literature.  We  recommend  simulating  the  Hugoniot  velocity  records 
measured  by  Larson  and  Anderson,  because  they  are  measured  histories,  not  derived  quanti¬ 
ties  like  the  peak  stress  values.  We  want  to  determine  if  good  correlations  with  the  velocity 


(1.2) 
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records  require  the  use  of  a  rate  effects  model,  and  if  so,  how  rate  effects  should  be  re¬ 
lated  to  pore  collapse.  If  the  rate  effects  and/or  stiffening  behavior  are  significant,  then  we 
recommend  performing  additional  SWAT  calculations  with  the  Mie-Gruniesen  model.  The 
appealing  aspect  of  the  Mie-Gruniesen  model  is  that  we  expect  it  to  simulate  higher  stress 
.  levels  near  the  SWAT  source  (high  pressure  regime)  than  the  rate-independent  constant  bulk 
modulus  model,  but  not  near  the  tunnels  (low  pressure  regime).  This  is  exactly  what  we 
accomplished  by  modeling  viscoplasticity  near  the  source  in  Figure  1-6,  but  not  near  the 
tunnels. 

Adding  brittle  damage  to  the  calculation  with  viscoplasticity  modeled  near  the  source  (but 
not  near  the  tunnels)  increases  the  late-time  springline  closures,  but  has  little  effect  on  the 
crown-invert  closures  or  the  velocity  and  stress  histories,  as  shown  in  Figure  1-8.  The  effect 
of  brittle  damage  on  the  springline  closures  is  most  noticeable  after  about  0.15  msec.  The 
brittle  damage  model  degrads  the  moduli  in  the  directions  of  principal  tensile  stress,  once 
an  initial  energy  threshold  is  exceeded.  Our  damage  model  is  loosely  based  on  research  by 
Ju  (1989).  Overall,  the  calculation  is  in  excellent  agreement  with  the  data. 

In  addition  to  the  damage  calculation  just  discussed,  we  ran  additional  brittle  damage  calcu¬ 
lations  with  and  without  modeling  viscoplasticity  near  the  source  and  tunnels,  but  the  trends 
were  always  similar.  Brittle  damage  either  increased,  or  had  little  effect,  on  the  late-time 
springline  closures.  In  some  cases,  reductions  in  the  late-time  stress  and  velocity  histories 
were  calculated  (after  about  0.10  msec). 

Adding  ductile  damage  to  the  brittle  damage  calculations  did  not  improve  comparisons  with 
the  measured  data,  with  or  without  modeling  viscoplasticity,  so  results  are  not  shown  here. 
For  example,  adding  ductile/brittle  damage  to  the  viscoplastic  calculation  previously  shown 
in  Figure  1-5  reduces  the  late-time  stress  histories,  in  poor  agreement  with  the  measured 
data,  and  increases  the  late-time  springline  closures.  The  ductile  damage  model  degrads  the 
moduli  in  the  directions  of  principal  compressive  stress,  once  an  initial  energy  threshold  is 
exceeded.  This  degradation  is  pressure-dependent,  so  substantial  damage  occurs  only  in  the 
low  pressure  regime  (below  about  25  MPa). 

The  reason  we  implemented  ductile  damage  in  addition  to  brittle  damage  is  that  our  brittle 
damage  implementation  did  not  simulate  a  "cross  effect”  between  degradation  in  one  princi¬ 
pal  direction  and  degradation  in  another,  as  suggested  by  Ortiz  (1985).  This  was  true  even 
when  we  attempted  to  model  brittle  damage  based  on  tensile  strain  directions  rather  than 
tensile  stress  directions®.  .An  example  of  modeling  the  cross  effect  is  an  unconfined  compres¬ 
sion  simulation  in  which  damage  initiates  in  the  lateral  directions  (due  to  the  presence  of 
lateral  tensile  strains),  causing  the  specimen  to  soften  in  the  axial  direction  (in  the  presence 
of  axial  compressive  strains).  VVe  were  unable  to  simulate  axial  softening  with  the  brittle 
model,  unless  ductile  damage  was  also  implemented.  Efforts  are  in  progress  to  develop  a 

®When  implementing  a  scalar  ani.sotropic  model,  one  must  choose  a  criteria  for  damage  evolution,  and 
define  the  projection  operators  (direction  of  damage).  For  example,  scalar  damage  criteria  are  often  based  on 
the  history  of  strains,  stress,  or  energy,  while  projection  operators  are  based  on  the  directions  of  maximum 
principal  strain  or  stress.  We  implemented  an  energy-based  damage  model  and  separately  examined  the  use 
of  projection  operators  based  on  maximum  principal  stress  and  maximum  principal  strain.  However,  neither 
implementation  simulated  a  "cross  effect.” 
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(a)  Closures  at  the  14.5  cm  tunnel. 
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(b)  Closures  at  the  19.2  cm  tunnel. 


Gages  PV2  PV5  PV8  PV11 


Time  (sec) 

(c)  Velocity  histories  at  four  gage  locations. 
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(d)  Stress  histories  at  two  gage  locations. 


Figure  1-8.  Modeling  brittle  damage  provides  more  late-time  springline  closure  than  the 
calculations  without  damage.  Here  we  modeled  viscoplasticity  near  the  source, 
but  not  near  the  tunnels. 
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(a)  Problem  definition.  (b)  Simulation  results. 

Figure  1-9.  This  unconfined  plane  strain  simulation  is  intended  to  represent  conditions  at 
the  tunnel  springline.  Note  that  little  softening  is  calculated  by  the  damage 
model. 


more  physically  based  anisotropic  model  which  includes  the  cross  effect. 

There  are  two  points  to  note  concerning  the  combined  effects  of  viscoplasticity  and  dam¬ 
age.  First,  viscoplasticity  increases  the  free-held  stress  histories,  while  damage,  if  anything, 
decreases  the  late-time  stress  histories.  Second,  modeling  viscoplasticity  near  the  tunnels 
tends  to  open  the  springlines,  while  modeling  damage  tends  to  close  the  late-time  spring¬ 
lines.  Perhaps  adjustments  in  our  damage  model,  or  our  hts  to  laboratory  data,  would  allow 
for  early-time  closure,  counteracting  the  opening  effects  of  modeling  viscoplasticity  near  the 
tunnels.  To  support  this  statement,  consider  the  single  element  plane  strain  simulation, 
unconhned  in  one  direction,  shown  in  Figure  1-9.  This  calculation  was  part  of  the  PTM 
calculators  benchmark  series  and  was  intended  to  simulate  the  conditions  at  the  springline. 
Note  that  our  model  predicts  very  little  softening,  particularly  as  the  strain-rate  increases. 
One  suggestion  is  that  modeling  more  severe  softening  might  result  in  early-time  springline 
closure,  in  better  agreement  with  the  measured  data. 

No  laboratory  test  data  is  availal)Ie  to  fit  the  softening  response  under  these  conditions.  It 
seems  reasonable  to  model  more  severe  softening  because  we  calculate  volume  expansion  in 
this  benchmark  problem.  In  general,  conventional  laboratory  tests  indicate  that  volume  ex¬ 
pansion  is  usually  accompanied  by  severe  softening.  This  benchmark  problem  and  additional 
SWAT  simulations  should  be  re-examined  once  a  more  physically  based  anisotropic  model  is 
developed. 

We  also  recommend  examining  and  displaying  the  calculated  contours  of  damage  near  the 
source  and  tunnels  for  comparison  with  the  post-test  swat  cross  sections  reported  b\’  Kiopp 
e:f.  ah  (199.5).  One  example  cross-section  is  reproduced  in  Figure  1-10.  SRI  concludes  that 
one  apparent  damage  mechani.sm  for  the  tunnels  is  shear  faults  and/or  cracks  extending  away 
from  the  springlines.  "Rubble  zones "  form  at  the  springlines  which  intrude  on,  the  tunnel  to 
cause  springline  closure.  We  interpret  these  crack  patterns  to  mean  that  modeling  damage 
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Figure  1-10.  SWAT  cross-section  showing  post-test  damage  near  the  source  and  tunnels. 

near  the  springlines  is  critical  to  obtaining  springline  closure.  Thorough  correlation  with 
the  SWAT  tests  should  include  measured  and  predicted  damage  patterns,  as  well  as  velocity, 
stress,  and  closure  histories. 

1.4  VISCOPLASTIC  MODEL  BEH.AVIOR. 

Modeling  rate  effects  with  viscoplasticity  significantly  affects  the  SWAT  computational  results, 
particularly  the  springline  closures.  Hence  it  is  important  to  understand  the  overall  behavior 
of  the  model  and  how  well  it  fits  available  test  data.  The  behavior  of  the  viscoplastic  model 
under  uniaxial  and  triaxial  loading  conditions  is  not  reported  in  the  literature,  so  a  overview 
is  presented  in  this  section. 

A  review  of  high  strain  rate  test  data  a\'ailable  in  the  literature  indicates  that  the  peak 
stress  obtained  from  Hopkinson  split  pressure  bar  tests  increases  with  strain  rate.  Data  for 
various  types  of  rock  from  Handin  and  IVie'dinan  (1976)  is  reproduced  in  Figure  1-11.  Note 
that  the  strength  for  Solenhofen  limestone  increa.ses  exponentially  with  log{c),  where  e  is  the 
strain  rate.  This  is  readily  apparent  above  rates  of  about  100/s.  On  the  other  hand,  some 
materials  like  basalt  and  tuff  exhibit  a  linear  relation  between  strength  and  log{e).  Hence  a 
rate-dependent  model  should  be  flexible  enough  to  capture  this  wide  range  in  behavior. 

Viscoplastic  models  are  typically  used  to  model  rate  effects  such  as  those  previously  shown 
in  Figure  1-11.  Duvaut-Lions  type  viscoplastic  models  are  easy  to  implement  into  existing 
cap  models.  In  general,  no  modifications  need  to  be  made  to  the  plasticity  algorithm;  the 
viscoplastic  implementation  is  separate  from,  and  follows,  the  plasticity  algorithm.  However, 
the  most  basic  Duvaut-Lions  model  does  not  provide  much  flexibility  in  fitting  data  because 
it  is  a  single  parameter  model.  The  original  constitutive  formulation  is  attributed  to  Duvaut 
and  Lions  (1972).  It  was  extended  to  three  dimensions  and  applied  to  multi-surface  cap 
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Figure  1-11.  High  strain  rate  data  for  rocks  indicating  an  increase  in  strength  with  strain 
rate  (Handin  and  Friedman,  1976). 

models  by  Simo  et.  al.  (1988).  Their  viscoplastic  model  requires  the  input  of  a  single 
parameter,  rj,  called  the  fluidity  parameter. 

1.4.1  Single-Parameter  Model. 

Some  e.xamples  of  good  and  poor  fits  to  laboratory  test  data  are  given  here  to  examine  the 
overall  behavior  of  the  single-parameter  Duvaut-Lions  model.  .A  good  fit  of  the  Duvaut- 
Lions  model  to  Solenhofen  limestone  data  from  Hopkinson  bar  tests  (Green.  1969)  is  shown 
in  Figure  1-12.  The  single  parameter  Duvaut-Lions  model  corresponds  to  n  =  1  solution: 
the  parameter  n  will  be  discussed  in  subsequent  paragraphs.  .Although  the  Duvaut-Lions 
model  fits  the  .Solenhofen  limestone  data  well,  it  would  not  provide  a  good  fit  to  all  rock 
data  previously  shown  in  Figure  1-11.  For  example,  a  poor  comparison  with  concrete  data 
is  shown  in  Figure  l-13b.  The  model  provides  more  strength  enhancement  in  tension  than 
in  compression,  in  qualitative  agreement  with  the  data.  However,  the  quantitative  fit  of  the 
model  to  data  is  poor  for  this  material.  The  Duvaut-Lions  curves  are  ’steeper'  than  the  data 
at  high  strain  rates.  The  concrete  data  hgure  was  reproduced  from  Ross  et.  al  (1992).  The 
data  of  interest  are  represented  by  the  symbols:  the  lines  represent  various  analyses  from 
(Ross.  1992)  that  should  be  ignored  for  our  purposes. 

.A  limited  amount  of  high  strain  rate  data  is  available  in  the  literature  for  stress  states  other 
than  uniaxial  stress  or  strain.  Solenhofen  limestone  data  (Brown.  1972)  for  a  Hopkinson 
bar  enclosed  in  a  pressure  vessel  are  reproduced  in  Figure  1-14.  There  is  a  55%  increase  in 
strength  at  a  conhning  pressure  of  30  .\IPa.  but  only  a  38%.  increase  in  strength  at  a  confining 
pressure  of  10  .MPa.  The  rate  effect  increases  with  confining  pressure  for  this  material. 
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Figure  1-12.  Good  fit  of  the  single  parameter  Duvaut-Lions  (n  =  1)  model  to  Hopkinson 
bar  data  for  Solenhofen  limestone. 

The  general  behavior  of  the  Duvaut-Lions  cap  model  implementation  under  triaxial  stress 
conditions  is  shown  in  Figure  1-15.  Without  the  cap,  the  peak  stress  predicted  with  the 
viscoplastic  model  appears  to  lie  on  a  'shifted'  plasticity  surface.  The  square  symbols  are  the 
peak  stresses  attained  in  viscous  TXr  simulations  with  confining  pressures  of  0,  25,  50,  and 
100  MPa.  These  same  computational  points  are  replotted  in  Figure  l-16a  versus  confining 
pressure  rather  than  pressure,  and  are  labeled  'dynamic’.  Also  plotted  for  comparison  are  the 
corresponding  strengths  obtained  with  rate-independent  simulations,  labeled  ‘static’.  This 
figure  indicates  that  as  the  confining  prcssun  increases,  the  predicted  rate  effect  decreases 
on  a  percentage  basis.  This  trend  is  more  readily  apparent  in  Figure  l-16b  in  which  the 
computational  points  are  plotted  as  the  ratio  of  the  dynamic  to  static  strengths.  This  is 
opposite  the  cpialitative  trend  previously  shown  for  Solenhofen  limestone  data  in  Figure  1- 
14. 

One  point  to  note  is  that  although  the  p(’ak  stress  appears  to  lie  on  a  shifted  plasticity  surface, 
the  stress-strain  curve  would  not  be  the  same  as  that  predicted  with  a  shifted  surface.  The 
stress  state  calculated  with  viscoplasticity  lies  outside  the  failure  surface.  The  stress-strain 
curve  calculated  with  viscoplasticity  is  nonlinear  once  the  stress  state  exceeds  the  plasticity 
surface  (which  does  not  shift),  whereas  the  stress-strain  curve  calculated  by  shifting  the 
plasticity  surface  would  be  linear. 
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Figure  I- 13.  Poor  comparison  of  the  single-parameter  Duvaut-Lions  smooth  cap 
implementation  to  concrete  data  (n  =  1). 
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Figure  1-14.  The  rate  effect  increases  with  confining  pressure  for  Solenhofen  limestone  data 
(Hopkinson  bar  enclosed  in  a  pressure  vessel). 


1.4.2  Strain-Rate  Dependent  Extension. 

These  comparisons  indicate  that  not  all  aspects  of  high  strain  rate  data  are  fit  well  and 
with  the  desired  flexibility.  To  provide  more  flexibility  in  fitting  the  data,  we  extended  the 
Du vaut- Lions  implementation  to  two  parameters,  effectively  making  the  fluidity  parameter 
a  function  of  strain  rate.  As  a  result  of  this  extension,  we  obtained  good  comparisons  with 
tensile  and  compressive  Hopkinson  bar  data,  as  shown  in  Figure  1-17.  .^n  overview  of  the 
extension  is  given  in  Appendix  B. 
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Figure  1-15.  The  peak  stress  attained  with  the  Duvaut-Lions  smooth  cap  implementation 
appears  to  lie  on  a  ’shifted'  plasticity  surface. 
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Figure  1*17.  Good  comparisons  of  the  strain-rate  dependent  viscoplastic  model  with 
concrete  data  in  =  5). 


(a)  Pressure-independent  model 


(b)  Pressure-dependent  model 


Figure  1-18.  The  pressure-dependent  viscoplastic  model  fits  Salem  limestone  Hugoniot  data 
better  than  the  pressure-independent  model 


1.4.3  Pressure-Dependent  Extension. 

In  general  the  original  and  modified  Duvaut-Lion  models  just  described  provide  a  good  fit  to 
Hopkinson  bar  data  for  a  variety  of  geological  materials,  but  not  to  Hugoniot  data  for  Salem 
limestone,  as  shown  in  Figure  1-lSa.  To  provide  a  better  fit,  we  added  pressure  dependence  to 
the  strain-rate  dependent  fluidity  parameter  to  model  greater  rate  effects  at  higher  pressures. 
Hence  we  are  effectively  modeling  greater  rate  effects  on  the  cap  than  on  the  shear  plasticity 
surface.  Good  agreement  between  the  pressure-dependent  model  and  Hugoniot  data  for 
Salem  limestone  is  shown  in  Figure  1-lSb,  Details  of  the  pressure-dependent  model  are 
given  by  Murray  and  Lindberg  (1996). 

One  might  think  that  we  could  model  different  rate  effects  on  the  cap  and  shear  failure 
surface  by  implementing  separate  pressure-independent  viscoplastic  formulations  for  the  cap 
and  shear  failure  surface.  Although  this  would  allow  one  to  separately  fit  the  Hugoniot  data 
and  Hopkinson  bar  data,  a  desirable  feature,  it  would  not  provide  a  better  qualitative  fit  to 
the  Hugoniot  data  shown  in  Figure  l-lSa.  A  pressure-independent  fluidity  parameter  would 
remain  constant  through-out  each  simulation.  Hence  we  could  obtain  a  different  fit  than 
previously  shown  in  Figure  1-lSa.  but  it  would  not  be  a  better  qualitative  fit.  To  obtain 
a  better  fit  to  the  Hugoniot  da.ta.  we  need  to  model  less  rate  effects  at  low  axial  strain 
(or  pressure)  and  more  rate  effects  at  high  axial  strain  (or  pressure)  than  modeled  with 
the  pressure-independent  formulation  in  Figure  1-iSa.  This  is  accomplished  by  varying  the 
fiuidst}*  parameter  as  a  function  of  pressure. 
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1.4.4  Inertial  Effects. 


One  current  issue  is  whether  the  sharp  increase  in  strength  with  strain  rate  commonly 
observed  in  Hopkinson  bar  tests  is  a  rate  effect  or  inertial  effect.  In  fact,  different  types 
of  high  strain  rate  tests  give  different  results  at  high  strain  rates,  indicating  that  inertial 
effects  may  be  less  significant  in  some  tests  than  others.  For  example.  Green  and  Perkins 
(1969)  used  a  compressive  split  Hopkinson  bar,  and  measured  a  100%  increase  in  strength  for 
Solenhofen  limestone  at  loading  rates  of  about  5000  s~^.  Lipkin,  Grady,  and  Cambell  (1977) 
used  a  torsional  split  Hopkinson  bar,  and  also  meeisured  about  a  100%  increase  in  shear 
strength  at  loading  rates  of  about  600  s~^.  On  the  other  hand.  Young  and  Dubugnon  (1977) 
used  a  reflected  shear  wave  technique,  and  determined  that  the  dynamic  shear  strength  for 
Solenhofen  limestone,  at  loading  rates  up  to  10®  s“^,  did  not  exceed  the  static  strength  by 
more  than  25%. 

Two  approaches  for  resolving  this  issue  are  computationally  (mesh  up  a  Hopkinson  bar  and 
predict  its  response)  and  experimentally.  Lateral  accelerations  of  a  Hopkinson  bar  specimen 
were  measured  by  Malvern  et.  al.  (1986)  and  Tang  et.  al.  (1992).  Malvern  concluded  that 
lateral  inertia  confinement  becomes  significant  at  strain  rates  of  about  100  s~^.  However, 
Tang  concluded  that  lateral  accelerations  were  too  small  to  account  for  the  factor  of  two 
enhancement  of  the  dynamic  compressive  strength  above  the  static  strength,  suggesting  that 
inertial  effects  are  insignificant  for  strain  rates  on  the  order  of  100/s.  This  is  an  issue  that 
should  be  resolved  in  future  efforts. 
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SECTION  2 

JOINT  MODEL  AND  STACK  OF  BRICK  SIMULATIONS 

2.1  JOINT  MODEL  OVERVIEW. 

We  use  a  rate-dependent  slideline  approach  for  explicitly  modeling  rock  joints.  To  provide 
a  better  fit  to  joint  normal  compressibility  data,  we  modified  the  DYNA3D  slideline  formu¬ 
lation  to  make  the  normal  stiffness  a  function  of  the  penetration  depth.  A  good  fit  of  the 
nonlinear  joint  model  to  normal  compressibility  data  is  shown  in  Figure  2-1.  To  provide 
more  flexibility  in  fitting  joint  shear  data,  we  modified  the  formulation  to  make  the  sheaj 
stiffness  independent  of  the  normal  stiffness.  The  joint  model  also  includes  a  rate-dependent 
Coulomb  friction  model  in  shear,  as  discussed  for  slidelines  in  the  DYNA3D  Users  Manual 
(Whirley  and  Engelmann,  1993) 

2.2  JOINT  RATE  EFFECTS. 

All  PTM  calculators  ran  a  benchmark  problem  to  exercise  the  joint  in  shear.  Two  blocks 
of  intact  rock  separated  by  a  diagonal  joint  are  subjected  to  normal  displacements.  The 
displacement  path  from  Simons  (1993)  is  reproduced  in  Figure  2-2.  Each  triangular  block 
was  modeled  with  a  single  element.  Computational  results  with  and  without  joint  rate 
effects  are  given  in  Figure  2-3.  They  indicate  that  the  joint  rate  model  lowers  the  joint 
shear  stress  at  high  strain  rates.  No  measured  data  was  available  for  verifying  our  joint  rate 
model.  Future  efforts  should  examine  available  test  data  to  determine  if  joint  rate  effects  are 
important,  and,  if  so,  determine  the  availability  of  joint  rate  data  for  fitting  and  validating 
the  model. 
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Figure  2-1.  The  nonlinear  joint  model  fits  normal  compressibility  data  better  than  the  linear 
model. 
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Figure  2-2.  Joint  shear  benchmark  geometry  and  loading. 
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Figure  2-3.  Preliminary  joint  sliding  calculations  with  and  without  rate  effects  (,/3  =  100 
indicates  rate  effects). 

2.3  ST.ACK  OF  BRICKS  COMPUTATIONAL  SU.VIMARV. 

The  Mighty  North  1  fMN'-l)  phenomenology  test  was  designed  to  investigate  the  effects  of 
joints  on  tunnel  deformation,  and  to  provide  data  for  benchmarking  and  model  validation.  It 
consists  of  a  tunnel  in  jointed  limestone  subjected  to  a  cylindrically-divergent  ground  shock 
load.  The  compressional  properties  and  damage  mode  of  the  jointed  limestone  "stack  of 
bricks”  (sOB)  arrangement  were  measured  by  AR.A.  for  benchmarking  calculations  (Chitty 
and  Blouin.  1995).  The  stack  of  bricks  test  configuration  is  given  in  Figure  2-4. 

Using  our  linear  joint  and  intact  rock  models,  we  performed  some  preliminary  plane  strain 
calculations  to  verify  our  SOB  generation  scheme  and  intact/joint  model  behaviors.  .Although 
the  cross-section  of  each  brick  is  nominallv  2-inch  bv  2-inch,  slight  variations  in  the  size  of 
each  brick  cause  vertical  and  horizontal  gaps  to  form  in  the  SOB.  Using  the  mesh  generator 
P.ATR.AX.  we  successfully  developed  a  system  to  scale  each  brick  model  to  form  a  stack 
of  irregular  bricks.  P.4TR.AN  conveniently  allows  the  user  to  create  one  brick,  then  scale 
and  translate  each  brick  to  form  a  stack  of  irregular  bricks.  The  scale  factors  are  based 
on  a  random  number  from  a  normal  distribution,  and  an  assumed  standard  deviation.  The 
translate  factors  are  calculated  by  placing  each  brick  on  the  highest  of  the  two  bricks  beneath 
it.  The  bottom  row  of  bricks  sits  on  a  steel  platten.  with  no  gaps  between  the  bricks  and 
platten. 


Load 


Figure  2-4.  Stack  of  bricks  layup  and  gage  locations. 

Our  preliminary  plane  strain  calculations  indicate  that  the  predicted  damage  pattern  for  the 
compressional  tests  is  sensitive  to  the  presence  of  irregular  gaps.  Damage  contours  calculated 
for  regular  (no  gaps)  and  irregular  (with  gaps)  bricks  are  shown  in  Figure  2-5  and  indicate 
that  modeling  irregular  gaps  causes  asymmetry  in  the  computed  damage  pattern.  Modeling 
irregular  gaps  also  produces  early-time  nonlinearities  in  the  stress-strain  response,  as  shown 
in  Figure  2-6.  Asymmetry  of  the  measured  damage  pattern,  and  early-time  nonlinearity  of 
the  measuied  stress-strain  response  were  also  observed  in  the  test. 
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(a)  Without  gaps. 


(b)  With  irregular  gaps. 


F'igurc  2-5.  hlodoliiig  irregular  gaps  produces  asAuumet rir  damage  pat  terns  in  the 
preliminary  "stack  of  l.u'icks**  simulation.  An  asymmetric  damage  pattern  was 
also  obser\'ed  in  the  stark  of  bricks  test. 


Figure  2-6.  The  preliminary  SOB  stress-strain  history  calculated  with  initial  gaps  is  more 
nonlinear  than  that  calculated  without  initial  gaps. 

These  plane  strain  calculations  were  intended  to  be  preliminary  runs  to  check  but  our  irreg¬ 
ular  SOB  generation  scheme  and  the  intact/joint  model  behaviors  before  running  more  time 
consuming  and  costly  three-dimensional  calculations.  To  complete  verification  of  our  joint 
model,  we  recommend  running  all  calculations  with  the  nonlinear  joint  model,  including 
horizontal  as  well  as  vertical  gaps,  and  analyzing  the  effect  of  boundary  conditions  on  the 
predicted  response,  including  rotation  of  the  plattens.  Another  possibility  is  setting  up  the 
irregular  brick  distribution  to  match  that  in  the  SOB  test.  ARA  measured  each  brick  dimen¬ 
sion  and  noted  its  placement  in  the  SOB.  In  this  way,  we  can  check  our  predicted  damage 
pattern  against  the  measured  damage  pattern  observed  in  the  test,  as  well  as  the  compare 
measured  versus  calculated  stress-strain  responses. 

Once  we  verify  our  model  for  this  particular  stack,  additional  parameter  studies  in  prepara¬ 
tion  for  MN-1  could  include: 


•  Three-dimensional  calculations  versus  plane  strain  calculations. 

•  The  effects  of  varying  joint  shear  stiffness,  as  well  as  varying  joint  normal  stiffness,  on 
the  predicted  response. 

•  The  effects  of  average  gap  size,  gap  standard  deviation,  and  brick  stacking  pattern  on 
stack  strength.  We  want  to  determine  if  stack  strength  is  primarily  determined  by 
average  gap  size  and  standard  deviation  regardless  of  the  stacking  pattern,  or  if  the 
stacking  pattern  plays  a  significant  role  in  strength.  This  is  important  because  we  will 
not  know  the  stacking  pattern  in  the  MN-1  test,  only  the  average  gap  size. 

•  The  effects  of  a  joint  rate  model  on  the  predicted  response  at  high  strain  rates. 
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SECTION  3 

SUMMARY  AND  LESSONS  LEARNED 


The  three-invariant  smooth-cap  model  with  damage  provides  good  fits  to  quasi-static  and 
high  strain  rate  data  for  limestone,  and  correlates  well  with  SWAT  data.  However,  the 
predictive  capability  of  the  model  has  not  been  established  for  tunnel  response  applications, 
because  only  a  few  static  and  dynamic  laboratory  tests  have  been  analyzed.  To  establish  a 
predictive  capability,  we  recommend  the  following  iterative  procedure: 

1.  Enhance  the  formulation  of  the  model.  At  this  time,  this  includes  additional  damage 
model  development,  understanding  the  differences  between  inertial  and  rate  effects, 
and  examining  alternatives  to  viscoplasticity. 

2.  Verify  the  formulation  through  correlations  with  dynamic  tunnel  test  data,  such  as 
SWAT.  We  recommend  comparing  contours  of  calculated  damage  with  post-test  obser¬ 
vations  of  the  fault/cracking  pattern,  as  well  as  comparing  measured  and  calculated 
free-field  and  tunnel  closure  histories.  Additional  parametric  studies  should  also  be 
performed  to  determine  the  sensitivity  of  tunnel  response  to  model  parameters  or  spe¬ 
cific  formulations  (like  damage  versus  rate  effects).  At  this  time,  we  reserve  judgement 
on  the  sensitivity  of  dynamic  tunnel  response  to  model  parameters  until  the  damage 
modeling,  rate  effects,  and  inertial  effects  issues  are  more  clearly  resolved. 

3.  Perform  pre-test  predictions  and/or  blind  post^test  predictions  of  dynamic  tunnel  tests 
performed  in  the  laboratory  or  field.  Accurate  predictions  are  the  only  means  of  vali¬ 
dating  a  code.  Post-test  correlations  should  only  be  used  for  identifying  improvements 
to  the  physical  basis  and/or  implementation  of  the  model. 

No  information  was  available  in  the  literature  describing  the  behavior  of  the  single-parameter 
Duvaut-Lions  viscoplastic  model  under  general  loading  conditions.  To  understand  the  model, 
we  examined  its  behavior  under  conditions  of  uniaxial  stress,  uniaxial  strain,  and  triaxial 
compression,  and  compared  the  behavior  to  geological/concrete  data  available  in  the  litera¬ 
ture. 

We  learned  that  the  single-parameter  Duvaut-Lions  model  (viscoplastic)  does  not  provide 
adequate  flexibility  in  fitting  a  wide  variety  of  high  strain  rate  geological  data  over  various 
regimes  (uniaxial  stress,  uniaxial  strain,  triaxial  stress).  We  improved  the  models  flexibility 
in  fitting  Hopkinson  bar  and  Hugoniot  data  by  extending  the  model  to  multiple  parameters. 
These  extensions  make  the  effective  fluidity  parameter  a  function  of  pressure  and  strain- 
rate.  The  extensions  are  based,  in  part,  on  analytical  comparisons  between  Duvaut-Lions 
and  Perzyna-type  models. 

Peak  velocities,  stresses,  and  closures  calculated  without  rate  effects  were  low  in  comparison 
with  measured  test  data.  To  improve  comparisons  with  the  test  data,  we  used  the  extended 
viscoplastic  model  in  our  SWAT  calculations.  Although  our  calculated  free-field  stress  and  ve¬ 
locity  histories  compared  well  with  the  measured  histories,  the  calculated  springlines  opened 
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slightly  whereas  the  measured  springlines  closed  slightly.  The  only  way  we  were  able  to  cal¬ 
culate  springline  closure  was  to  model  rate  effects  near  the  source,  but  not  near  the  tunnels. 
We  note  that,  near  the  source,  the  pressure  is  high  and  the  stress  state  lies  on  the  cap.  Near 
the  tunnels,  the  pressure  is  low  and  the  stress  path  encounters  the  shear  failure  surface. 

Two  approaches  could  be  taken  to  attempt  better  agreement  with  the  springline  closures. 
The  first  approach  is  to  continue  to  tune  our  fits  and/or  viscoplastic  formulation.  The 
general  approach  of  the  PTM  group  calculators  was  to  model  different  rate  effects  on  the  cap 
than  on  the  shear  failure  surface,  z.e.  strong  rate  effects  on  the  cap,  but  minimal  rate  effects 
on  the  shear  failure  surface.  The  assumption  here  is  that  different  physical  processes  are 
taking  place:  volume  crushing  on  the  cap,  and  shear  failure  on  the  failure  surface.  However, 
two  modifications  were  made  to  the  viscoplastic  model  to  get  it  to  agree  with  laboratory 
data  suggesting  that  viscoplasticity  may  not  be  the  best  method  of  modeling  these  different 
mechanisms.  In  fact,  the  general  behavior  of  the  original  viscoplastic  (before  modifications) 
was  that  rate  effects  decreased  with  pressure,  at  odds  with  the  desired  trend. 

The  second  approach  is  to  examine  alternatives  to  viscoplasticity.  One  example  is  use  a 
nonlinear  Mie-Gruniesen  model  in  place  of  a  constant  bulk  modulus  model.  Comparisons 
with  measured  particle  velocity  records  from  Hugoniot  tests  should  indicate  whether  pore 
collapse  rate  effects  still  need  to  be  modeled.  Although  no  data  was  available  for  Salem  lime¬ 
stone  at  very  high  compaction  levels,  many  geological/concrete  materials  exhibit  a  stiffening 
behavior  at  high  pressure.  Such  behavior  is  simulated  by  the  Mie-Gruniesen  model,  but 
not  by  the  constant  bulk  modulus  model.  We  recommend  this  approach  for  future  tunnel 
response  calculations,  at  least  for  comparison  with  viscoplastic  model  results. 

These  efforts  show  that  no  single  set  of  laboratory  test  data  can  be  used  to  formulate  and 
verify  a  comprehensive,  three-dimensional  rate-dependent  model.  We  want  to  formulate 
models  with  the  right  trends  over  a  range  in  loading  regimes.  This  can  be  accomplished 
by  examining  a  wide  variety  of  geological/concrete  materials  and  tests,  such  as  tensile  and 
compressive  Hopkinson  bar,  plate-slap  Hugoniot  tests,  and  divergent  flow  tests.  Once  the 
general  behavior  of  a  model  is  unchnstood.  then  the  specific  behavior  of  the  model  can  be 
verified  for  a  material  of  interest,  such  as  Salem  limestone.  This  same  approach,  general 
formulation  then  specific  fits,  can  \)r  used  wlien  developing  damage  models. 

We  implemented  a  scalar  anisot  ropi(*  damage'  mudt'l  in  order  to  to  simulate  strain-softening 
and  modulus  reduction  in  th('  tensih*  and  low-confining  pressure  regimes.  The  model  is 
loosely  based  on  research  by  ,]u  The  motivation  for  implementing  an  anisotropic 

model  instead  of  an  isotroi)ic^  model  is  that  damage  in  geological  materials  is  directional 
in  nature.  Anisotropy  adds  more'  i)hysics  to  the  model,  without  requiring  any  additional 
laboratory  data  to  fit  the  model. 

We  used  the  anisotropic  damage  model  in  our  SWAT  calculations,  with  and  without  modeling 
viscoplasticity.  We  obtained  our  best  correlations  with  the  measured  SWAT  histories  using 

^Interested  readers  should  refer  to  Murray  and  Lewis  (1995)  for  a  description  of  the  isotropic  damage 
formulation  that  aptek  previously  implemented  into  dyna3d  on  dna’s  Conventional  Weapons  Effects  (sc 
ewe)  Program.  This  model  ha.s  sucessfully  been  used  to  predict  damage  to  reinforced  concrete  bunkers 
subjected  to  internal  detonations. 
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a  combination  of  brittle  damage  and  viscoplasticity  modeled  near  the  source.  Modeling 
damage  and  viscoplasticity  tend  to  have  counteracting  effects  on  some  of  the  calculated 
histories.  In  general,  viscoplasticity  increases  the  free-held  stress  histories,  while  damage 
decreases  the  late-time  stress  histories.  Modeling  viscoplasticity  near  the  tunnels  opens  the 
springlines,  while  modeling  damage  tends  to  close  the  springlines,  at  least  at  late-time. 

Although  our  SWAT  correlations  are  reasonable,  we  recommend  further  anisotropic  damage 
model  development  for  two  reasons.  First,  the  model  simulates  little  strain-softening  under 
conditions  of  “unconhned  plane  strain,”  such  as  occurs  near  the  springline  closure.  We 
would  like  to  determine  if  more  severe  softening  under  these  conditions  allows  for  early- 
time,  rather  than  late-time,  springline  closure,  in  better  agreement  with  the  opening-closing 
response  of  the  measured  data.  Second,  the  model  does  not  simulate  a  “cross  effect’  between 
degradation  in  one  principal  direction  and  degradation  in  another,  hence  the  physical  basis 
of  our  implementation  is  not  satisfactory.  Future  efforts  should  examine  the  combinations 
of  scalar  damage  criteria  and  projection  operators  that  will  realistically  produce  a  cross 
effect.  Once  the  general  anisotropic  formulation  is  verified  against  test  data  available  in  the 
literature,  additional  SWAT  parametric  studies  should  be  performed,  with  and  without  the 
proposed  Hugoniot  and  viscoplastic  models. 

At  the  beginning  of  the  current  program,  predictions  were  made  with  geological  models  with 
limited  physics;  damage  (softening  and  modulus  reduction)  was  neglected  while  rate  effects 
were  modeled  in  a  very  approximate  manner.  Thus  most  differences  between  the  measured 
and  calculated  histories  were  initially  attributed  to  how  we  modeled  rate  effects.  The  result  of 
this  sequential  approach  is  that  modeling  rate  effects  received  more  attention  than  modeling 
damage  effects.  However,  subsequent  damage  model  calculations  and  post-test  observations 
of  the  SWAT  fault/cracking  patterns,  indicate  that  modeling  damage  may  play  an  important 
role  in  correlating  with  the  SWAT  data.  For  future  efforts  we  recommend  a  more  global 
approach  in  which  known  physical  mechanisms  are  modeled  before  parametric  studies  are 
performed  and  model  sensitivities  are  concluded. 

Kfforts  to  develop  and  verify  a  jointed  rock  model  are  still  in  progress.  We  sucessfully 
developed  a  nonlinear  joint  model  which  fits  normal  compressibility  data  well.  Preliminary 
plane-strain  "stack  of  brick"  simulations  were  |)erformed  to  understand  the  effects  of  joint 
gaps  on  predicted  damage  patterns  and  stress-strain  curves.  Irregular  gaps  cause  asymmetry 
in  the  predicted  damage  pattern  and  nonliiu'arit y  of  the  stress-strain  response.  Future  efforts 
should  include  three-dimensional  calculations  for  comparison  with  the  measured  stress-strain 
curves  and  post-test  damage  pattern,  to  verify  our  jointed  rock  model.  However,  as  previously 
stated,  accurate  pi'fdictions  are  the  only  iiu'ans  of  validating  a  code.  Future  calculations 
should  include  “blind"  post-test  .\IN-1  predictions  and  field  test  predictions.  Unlike  the 
manufactured  joints  in  the  MN-1  tests,  natural  rock  joint  patterns  will  be  random  (and 
unknown)  in  field  tests,  which  will  be  a  challenge  for  future  efforts. 
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APPENDIX  A 


FITS  OF  THE  INTACT  ROCK  MODEL  TO  DATA 


The  three-invariant  smooth-cap  model  with  brittle  and  ductile  damage  provides  good  fits  to 
quasi-static  laboratory  test  data  for  Salem  limestone  as  demonstrated  in  Figures  A-1  and 
A-2.  The  data  were  provided  by  ARA  (Chitty  and  Blouin,  1995)  and  RE/spec  (Fossum 
et.  aL,  1995).  Both  the  peak  stress  and  post-peak  softening  responses  fit  the  data  well. 

No  high  strain  rate  laboratory  tests  were  conducted  on  Salem  limestone  by  ARA  or  RE/spEC, 
although  Hugoniot  data  for  Salem  limestone  was  available  in  the  literature  (Larson,  1979). 
The  fit  of  our  pressure-dependent  viscoplastic  model  to  this  data  is  shown  in  Figure  A-3a. 
As  an  alternate  specification,  we  also  fit  our  viscoplastic  model  directly  to  SWAT  data,  as 
discussed  in  Section  1.3.  Stress-strain  curves  resulting  from  this  fit  are  shown  in  Figure  A-3b. 
They  indicate  greater  rate  effects  than  the  fit  to  Hugoniot  data  shown  in  Figure  A-3a. 

The  SWAT  and  Hugoniot  data  is  primarily  high  pressure  data,  so  the  rate-dependent  behavior 
in  uniaxial  stress  still  has  to  be  specified.  Based  on  a  general  review  of  Hopkinson  bar  data 
available  in  the  literature,  we  assumed  that  strength  enhancement  increases  by  about  a 
factor  of  two  at  high  strain  rates  near  lO"*.  This  is  true  for  Solenhofen  limestone  shown  in 
Figure  1-12  and  reported  by  Green  (1969)  These  assumed  ‘fits’  are  given  in  Figures  A-4  and 
A-5. 
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Figure  A-1.  Good  fits  of  model  to  triaxial  extension  data  for  Salem  limestone. 
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Figure  A-.3.  Uniaxial  strain  behavior  of  the  model  at  four  strain  rates  for  two  different  fits 
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(a)  Unconfined  compression. 
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Figure  A-4.  Assumed  strength  versus  strain  rate  behavior  of  model  showing  greater  rate 
effects  in  tension  than  in  compression. 
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Figure  A-5.  Uniaxial  stress  behavior  of  the  model  at  four  strain  rates. 


APPENDIX  B 


ANALYTICAL  SOLUTIONS  FOR  DUVAUT-LION  AND  PERZYNA  MODELS 


We  use  the  viscoplastic  update  scheme  suggested  by  Simo  et.  al.  (1988),  except  we  make 
the  effective  fluidity  parameter,  a  function  of  strain  rate,  as  follows; 


/^\ 

r}n  =  Vi  jj 
1 

At/r]n  +  l 

i,,  =  (1  -  a)iij  +  adfj 


(B.l) 


(B.2) 


rjn  —>  0  (T,j  =:  (T,j  Inviscid  Solution 
rjn  oo  aij  —  Elastic  Solution 

Here  Oij  and  a,j  are  the  viscid  and  inviscid  stress  tensors,  respectively,  and  dj”®'  is  the  elas¬ 
tic  trial  stress.  At  each  time  step,  and  are  calculated  and  saved  by  the  plasticity 

algorithm.  Then  we  apply  Equation  B.2  to  update  the  viscous  solution.  Note  that  the  invis¬ 
cid  and  instantaneous  elastic  solutions  arc  obtained  as  limiting  cases.  The  single  parameter 
Duvaut-Lions  model  corresponds  to  the  ??  =  1  solution. 

One  motive  for  extending  the  single-parameter  Duvaut-Lion  model  to  two  parameters  is  to 
provide  a  more  flexibility  in  fitting  high  strain  rate  data.  Similarly,  Simo  et.  al.  (1988)  suggest 
a  form  for  exponential  and  power  law  models  which  are  more  flexible  in  fitting  data  than  their 
basic  single  parameter  model.  However,  these  more  elaborate  models  require  the  solution  of  a 
nonlinear  equation  at  each  and  every  time  step  of  the  finite  element  calculation.  The  solution 
of  the  nonlinear  equation  provides  the  interpolation  parameter  a  as  an  implicit  function  of 
strain  rate.  On  the  other  hand,  our  proposed  implementation  provides  the  interpolation 
parameter  o  as  an  explicit  function  of  strain  rate.  Our  implementation  eliminates  the  need 
to  solve  a  nonlinear  equation  at  each  time  step.  This  is  an  advantage  when  running  large 
calculations  where  CPU-time  is  an  issin'.  Ih'iice.  a  second  motive  for  extending  the  single¬ 
parameter  Duvaut-Lion  model  to  two  i)arametcrs  is  to  develop  an  efficient  run-time  model. 

Here  we’ve  converted  the  single-parameter  Duvaut-Lion  model  to  a  two-parameter  model. 
We  say  two  parameters,  rather  than  three  parameters,  because  t/i  can  be  regrouped  with  7 
to  yield  one  less  parameter.  The  two  parameters  for  fitting  Perzyna-type  behavior  are  n  and 
7  with  7i  =  1.  The  one  parameter  for  fitting  Duvaut-Lions  type  behavior  is  71  with  n  =  1. 

The  form  of  the  fluidity  parameter  in  Equation  B.l  is  based  on  comparisons  between  ana¬ 
lytical  solutions  for  Duvaut-Lions  and  Perzyna  power  law  models.  A  thorough  discussion  of 
the  derivation  is  given  by  Murray  and  Lindberg  (1996).  Perzyna-type  power  law  models  are 
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more  flexible  in  fitting  data  than  the  Duvaut-Lions  model,  because  they  are  three-parameter 
models  rather  than  single  parameter  models.  However,  they  are  not  as  easy  to  implement 
into  existing  cap  models.  In  fact,  the  plasticity  algorithm  must  be  rewritten  to  accommodate 
the  viscoplastic  implementation.  This  is  a  time-consuming  and  tedious  task,  and  hence  was 
not  undertaken.  See  Perzyna  (1971)  for  a  discussion  of  Perzyna  type  models. 

For  a  state  of  uniaxial  stress,  the  analytical  Duvaut-Lions  and  Perzyna  solutions  provide  the 
peak  dynamic  strength  (viscid  stress)  as  a  function  of  static  strength  (inviscid  stress)  and 
strain  rate  as  follows: 

viscid  —  ^inviscid  “I”  )  (^'^) 

where  E  is  Young’s  modulus.  One  can  use  Equation  B.3  to  estimate  the  effective  fluidity 
parameter  when  fitting  high  strain  rate  data.  The  resulting  behavior  in  Figure  B-1.  Here 
we’ve  plotted  the  viscid/inviscid  stress  ratio  as  a  function  of  e/7  for  three  values  of  n.  The 
n  —  1  solution  provides  the  largest  increase  in  strength  over  the  smallest  range  in  strain  rate: 
it  is  the  ‘steepest’  solution. 

These  numerical  peak  stress  solutions  may  be  fit  to  Hopkinson  bar  data  by  overlaying  the 
solutions  onto  data  plots  and  shifting  right  or  left  until  one  of  the  numerical  solutions  fits  the 
data,  as  previously  demonstrated  for  Solenhofen  limestone  data  (Green,  1969)  in  Figure  1-12. 
The  value  of  n  is  obtained  from  the  chosen  numerical  solution  and  the  value  of  7  (Perzyna 
Parameter)  is  obtained  by  equating  the  overlay  abscissa  (e/7)  to  the  data  abscissa  (e).  As 
previouslv  shown  in  Figure  1-12,  the  n  =  1  Duvaut-Lions  solution  provides  the  best  fit  to 
the  limestone  data,  although  the  lone  data  point  at  e  =  0.1  s"^  is  not  fit  well.  However,  the 
n  =  5  solution  provides  a  better  fit  to  the  concrete  data  than  the  n  =  1  solution  previously 
shown  in  Figure  1-13.  This  good  fit  was  previously  shown  in  Figure  1-17.  Our  two-parameter 
formulation  provides  Perzyna-like  flexibility  in  fitting  data  with  the  ease  of  implementation 
of  a  Duvaut-Lions  model. 
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(a)  Linear  abscissa. 


LOG(Strain  Rate  /  Perzyna  Parameter) 

(b)  Logarithmic  abscissa. 


Figure  B-1.  Numerically  calculated  strength  enhancement  indicating  that  the  n  =  1  solution 
has  the  steepest  slope. 
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